Für die meisten GIS Anwendungen in R könnt ihr das package “sf” benutezen. Eine ausführliche Anleitung zu fast allem, was sf zu bieten hat findet ihr hier. Im Rahmen des Projektes werdet ihr sf brauchen, um die einzelnen Probelstellen ihren Catchments zuzuweisen. In diesem Skript werden wir die Probestellen Bundesländern zuordnen, und uns die Catchments auf Detuschland beschränken.

Als erses laden wir das Package.

library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

Dann zuerst die Probendaten. Die Probedaten liegen im .RDS Format vor. Dieses Format ist speziell für R und alle R Objkete lassen schnell darin speichern saveRDS() und landen readRDS(). Besonders bei rumlichen Daten kann die größe der Dateien und die Ladedauer damit stark reduziert werden. Die Daten können allerdings von keinem anderen Programm (z.B. QGIS) geöffnet werden.

mzb <- readRDS("data/germany.rds")

und die Catchments. Die Catchments liegen als Shapefile vor. Shapefiles kennt ihr hoffentlich gut aus dem GIS Kurs. Sie bestehen immer aus mehreren Datein aber nur die .shp Datei kann genutzt werden um die räumlichen Daten zu öffnen. Um die Shapefile in R einzulesen nutzen wir die erste Funktion aus sf: sf_read. Alle Funktionen in sf fangen mit dem beiden Buchstaben st (kurz für Spatio-temporal) an. Da R Studio über eine autocomplete Funktion verfügt kann man also immer sf_ eintippen und bekommt die vollständige Liste an Funktionen aus dem sf Package.

cat <- st_read("data/lemm/MultipleStress_RiverEcoStatus.shp")
## Reading layer `MultipleStress_RiverEcoStatus' from data source 
##   `C:\Users\jonat\Documents\Uni\teaching\projekt_uwi\21_betadiversity\data\lemm\MultipleStress_RiverEcoStatus.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 52847 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -10.4779 ymin: 34.64526 xmax: 32.93388 ymax: 70.1431
## Geodetic CRS:  WGS 84

Einen Polygonlayer der deutschen Bundesländer können wir mit den geodata package runterladen. Den Code werde ich hier ich explizit erklären, mehr Infos zu dem geodate package findet ihr hier.

library(geodata)
## Lade nötiges Paket: terra
## terra version 1.4.20
ger1 <- gadm(country = "DEU", level = 1, path = "data")
ger1 <- st_as_sf(ger1)

Häufig benutzte Funktinen

Im folgenden Abschnitt werde ich einge Funktionen demonstrieren die man allgemein häufig gebrauchen kann. Die ersten Funktionen (filter(), select(), ) stammen aus dem dplyr package. Man benutzt sie um bestimmte Reihen (filter) oder Spalten (select) einer Tabelle auszuwählen. Die Probedaten habem viele Spalten (n = ncol(mzb)), die meisten davon sind für uns hier uninterissant. Wir wählen nur die site_id, den Genusnamen und das Datum aus. Wenn man select auf ein sf Objekt anwendet wird die Spalten mit den Koordinaten automatisch mit ausgewählt.

library(dplyr)
## 
## Attache Paket: 'dplyr'
## Die folgenden Objekte sind maskiert von 'package:terra':
## 
##     intersect, src, union
## Die folgenden Objekte sind maskiert von 'package:stats':
## 
##     filter, lag
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     intersect, setdiff, setequal, union
mzb2 <- select(mzb, c("site_id", "genus", "date"))
head(mzb2)
##   site_id     genus       date
## 1   00001   Antocha 2014-04-15
## 2   00001   Ancylus 2014-04-15
## 3   00001    Baetis 2014-04-15
## 4   00001      <NA> 2014-04-15
## 5   00001 Dicranota 2014-04-15
## 6   00001   Dugesia 2014-04-15

Aus diesem Datensatz können wir jetzt noch die Reihen auswählen die Baetis Arten enthalten. Wie immer in R müssen wir zwei Gleichzeichen benutzen um ist gleich auszudrücken, da ein einfaches Gleichzeichen auch für Objektzuweisungen genutzt werden kann.

mzb3 <- mzb2 |>  filter(genus == "Baetis")

Anbei ein weiteres beispiel in dem wir drei Bundsländer aus dem ger1 Objekt herrausfiltern. Das %in% Zeichen bedeuted ist Element von. Der Befehl sucht also alle Reihen aus ger1 aus in dennen die Spalte NAME_1 einem der Elemente des Vektors `c(“Bayern”, “Berlin”, “Bremen”)``entspricht.

ger2 <- filter(ger1, NAME_1 %in% c("Bayern", "Berlin", "Bremen"))

Das können wir uns auch auf einer Karte ansehen:

library(mapview)
mapview(ger2)

Manchmal möchte man auch nur eine Spalte entfernen oder alle Reihen auswählen die einen bestimmten Eintrag nicht haben. Das geht in R mit dem Ausrufezeichen. Beispiel:

cat2 <- select(cat, !eco_stat_n)
cat2 <- filter(cat, eco_stat_2 != "bad")
cat2 <- filter(cat, !eco_stat_2 %in% c("bad", "moderate"))

Probedaten in sf überführen.

Die Probedaten sind aktuell kein räumliches Objekt. Sie sind eine einfacher data.frame (eventuell wird bei euch data.table angezeigt). Um aber Probedaten auszuwählen, die in einem bestimmten Bundesland oder Catchment liegen müssen sie räumlich werden - also sf Objekte. Das machen wir mit der Funktion st_as_sf(). Dabei gibt es drei Sachen zu beachten. 1. Aktuell entspricht jede Reihe in mzb der Beobachtung von einer Art, an einem Datum und an einer Probestelle. Würden wir diese Tabelle so in eine sf Objekt umwandeln, hätten wir für jede Art and jeden Datum einen Punkt. Da viele der Punkte übereinander liegen würden bringt uns das nichts. Wir können den Rechenaufwand also reduzieren, indem wir nur eine Reihe pro Probestelle behalten. Auch dafür können wir eine Funktion aus dplyr nutzen. Mit distinct() wählen wir eine Reihe pro Eintrag der angegebenen Variable (site_id) aus.

mzb2 <- distinct_(mzb, "site_id", .keep_all = TRUE)
## Warning: `distinct_()` was deprecated in dplyr 0.7.0.
## Please use `distinct()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
  1. Wir müssen der Funktion sagen in welchen Spalten sie die Koordinaten findet. In diesem Fall x.coord und y.coord.
  2. Wir mpssen der Funktion sagen in welchem Koordinatenbezugssystem (KBS, crs in der Funktion) die Daten projeziert sind. Hier ist der EPSG code als Variable im Datensatz enthalten. Wir wählen nur den ersten Eintrag von mzb$EPSG aus da die Funktion verwirtt ist wenn wir einen Vektor als Argument für das Koordinatenbezugssystem angeben.
sf_mzb <- st_as_sf(mzb2, coords = c("x.coord", "y.coord"), crs = mzb$EPSG[1])

Letztlich wollen wir hier, als Beispiel, nur die Probestellen raussuchen, die in Bayern liegen. Wir müssen darauf achten, dass beide Objekte (die Probestellen und Bayern) im gleichen KBS liegen. Das können wir mit der Funktion st_crs prüfen.

st_crs(sf_mzb) == st_crs(ger1)
## [1] FALSE

Aktuell haben sie noch unterschiedliche KBS. Wir müssen also einen von beiden transformieren.

ger2 <- st_transform(ger1, crs = st_crs(sf_mzb))

Die Auswahl der Probestellen erfolgt durch st_filter(). Die Funktion filtert alle Objekte aus dem ersten Argument die innerhalb des zweiten Argumentes liegen. Ganz ähnlich funktioniert die Funktion st_join() mit der ihr Variablen von einem Layer (den Catchments) auf einen anderen (die Probestellen) übertragen könnt.

bav <- filter(ger2, NAME_1 == "Bayern" )
sf_mzb_bav <- st_filter(sf_mzb, bav)
mapview(sf_mzb_bav)

Für die Analysen würde es Sinn machen sich auf die aktivsten Jahre (1986:2013) und die Monate April bis Juni zu konzentrieren.

Um aus einem Datum, wie der date Spalte im Probedatensatz den Monat zu extrahieren könnt ihr die Funktion month() aus dem lubridate package verwenden.

library(lubridate)
## 
## Attache Paket: 'lubridate'
## Die folgenden Objekte sind maskiert von 'package:terra':
## 
##     intersect, union
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     date, intersect, setdiff, union
months <- month(sf_mzb_bav$date)

Die Catchments

Die Polygone der Shapefiles sind leider nicht gut gemacht. Was genau das Problem ist weis ich nicht, es könnte sein, das einige Polygone nicht perfekt schließen oder sich die Linien eines Polygons selber kreuzen. Wie dem auch sein all das können wir mit der wunderbaren Funktion st_make_valid() beheben.

cat2 <- st_make_valid(cat)

Wie oben erwähnt wollen wir nur die einzelnen Catchments die aktuell noch große Teile Europas bedecken auf Detuschland reduzieren. Das geht mit der bereits benutzen st_filter() Function.

cat2 <- st_transform(cat2, crs = st_crs(ger1)) 
cat3 <- st_filter(cat2, ger1)
mapview(cat3)

betapart

Für die Analyse der betadiversität müssen wir die so ausgewählten Daten allerdings noch vorbereiten.

Die Daten müssen in dem Format eine Reihe pro Probenahme und eine Spalte für jedes Taxon vorliegen. Auch hierfür müssen wir wieder ein neues Package verwenden: tidyr. Zuerst müssen wir allerdings die Spalte mit den Koordinaten, die immer mitgenommen wird auch wenn sie nicht explizit gewählt wird das geht mir st_drop_geometry.

library(tidyr)
## 
## Attache Paket: 'tidyr'
## Das folgende Objekt ist maskiert 'package:terra':
## 
##     extract
sf_mzb_bav <- st_drop_geometry(sf_mzb_bav)
focal_bio <- pivot_wider(sf_mzb_bav, id_cols = gr_sample_id, names_from = "lowest.taxon", values_from = abundance, values_fill = 0)
focal_bio <- select(focal_bio, !gr_sample_id)

pivot_wider() nimmt die Tabelle und formt sie um. Eine Reihe pro id_cols, eine Spalt pro names_from. Die Werte kommen aus der Spalte values_from und in Fällen in dennen die Kombination nicht existiert wird eine 0 (values_fill) eingesetzt.

Wir nutzen nur Presence-Absence und keine Abundanzen. Da die Daten aktuell noch Abundanzen sind müssen wir sie Umwandeln.

library(vegan)
## Lade nötiges Paket: permute
## Lade nötiges Paket: lattice
## This is vegan 2.5-7
focal_bio <- decostand(focal_bio, method = "pa")

Mit diesem Datensatz können wir jetzt Betadiversitäten berechnen.

library(betapart)
bc  <- betapart.core(focal_bio)
bc2 <- beta.multi(bc)

In bc2 ist nun die Betadiversität von Bayern … sozusagen. Es macht allerdings nicht wirklich Sinn das auf einer solchen Skala zu berechnen, das ist hier ja nur zur Demonstration.

Eine letzte Sache noch. Ihr müsste den Code nicht für jedes Catchment einzeln schreiben. Um euch Zeit zu sparen und euren Verstand zu bewaren könnt ihr einen for Loop benutzen. For Loops haben eine Variable (i) für die jede Runde ein anderer Wert angenommen wird. Mit dieser Variable wird dann der Code im Loop ausgeführt. Am Besten seht ihr das an einigen Beispielen:

# Der Code (print(i) wird für jeden Wert zwischen 1:10 evaluiert)
for (i in 1:10) {
        print(i)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
for (i in 1:10){
        x <- 2 + i
        print(x)
}
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
# Hier speicher ich die Ergebnisse jeder Runde
save_list <- list()

# Die Variable muss nicht i heißen ihr könnt sie nennen wie ihr wollt. 

for (unpraktischerlangername in c(5,10,15,20,30,40,50)){
        x <- filter(sf_mzb_bav, gr_sample_id %in% unique(sf_mzb_bav$gr_sample_id)[1:unpraktischerlangername])
        x <- pivot_wider(x, id_cols = gr_sample_id, names_from = "lowest.taxon", values_from = abundance, values_fill = 0)
        x <- select(x, !gr_sample_id)
        x <- decostand(x, method = "pa")
        bc  <- betapart.core(x)
        bc2 <- beta.multi(x)
        bc2 <- data.frame(SIM = bc2[[1]], 
                          SNE = bc2[[2]], 
                          SOR = bc2[[3]])
        save_list[[length(save_list) + 1]] <- bc2
}

x <-  do.call("rbind", save_list)